Genome and life-history evolution link bird diversification to the end-Cretaceous mass extinction

Complex patterns of genome evolution associated with the end-Cretaceous [Cretaceous-Paleogene (K–Pg)] mass extinction limit our understanding of the early evolutionary history of modern birds. Here, we analyzed patterns of avian molecular evolution and identified distinct macroevolutionary regimes across exons, introns, untranslated regions, and mitochondrial genomes. Bird clades originating near the K–Pg boundary exhibited numerous shifts in the mode of molecular evolution, suggesting a burst of genomic heterogeneity at this point in Earth’s history. These inferred shifts in substitution patterns were closely related to evolutionary shifts in developmental mode, adult body mass, and patterns of metabolic scaling. Our results suggest that the end-Cretaceous mass extinction triggered integrated patterns of evolution across avian genomes, physiology, and life history near the dawn of the modern bird radiation.


Supplementary Text
Molecular and life-history datasets Our final combined nuclear genetic dataset is a matrix of 400 rows (two alleles per taxon for each locus) and 682,178 aligned base pairs.These data reflect an effective increase in called bases by a factor of ~3.45 and an increase in alignment length by a factor of 1.73 compared with the original dataset from Prum et al. 2015 (37), accompanied by increases in the percentage of undetermined and gap characters of less than 10% overall.Only ~14% of the characters in the combined dataset are gaps "-", and ~8% are undetermined characters.The final combined mitochondrial dataset encompasses all 15 targeted regions (13 proteins and 2 rRNAs) and 200 taxa with an average alignment length of 13,769 base pairs (Table S3).Both life-history datasets reflect the phylogenetic sample of avian lineages in Prum et al. 2015 (37).For the first dataset, we assembled data for body mass (0% missing), modeled generation length (0% missing), latitude centroid (0% missing), mean clutch size (~13.6%missing), annual adult survival (77% missing), age at first breeding (58% missing), maximum known longevity (54.5% missing), and developmental mode ("ChickPC1" -0% missing).Categorical variables have no missing data (Fig. 2, left).For the second dataset highlighting avian metabolism, our procedure to identify phylogenetically equivalent swaps generated ~70% coverage for BMR and 100% for body mass (Table S4a).

Molecular dataset details
The percentage of undetermined 'N' characters was relatively low in exons (~1.3%) but substantially higher in introns (~24%) and UTRs (~18%).Our final filtered dataset comprised 384 exons (average ~1,198bp), 365 introns (average ~404bp), 147 UTRs (average ~453bp), and 14 loci which were noncoding but did not clearly map to a non-coding genic region (average ~578bp).The reprocessed dataset also reflects an increase in the number of parsimony informative sites by a factor of 2.18 for the phased alignment (~1.8 for a single allele) relative to the Prum et al. 2015 (37) dataset.These encouraging results imply that many previously generated phylogenomics datasets may be effectively reprocessed with updated approaches to generate significantly larger and well-annotated resources for studying phylogeny and molecular evolution.Due to computational bottlenecks, our molecular model shift analyses were conducted on sub-sampled sequence datasets representing one sampled allele from each species (we provide the re-assembled datasets as supplementary data).Lastly, the final combined mitochondrial dataset encompassed all 15 targeted regions (13 proteins and 2 rRNAs) and covered 200 taxa with an alignment length of 13,769 base pairs (Table S3).These data also have a low percentage of missing or undetermined characters (12.8%).
Assessing the coincidence of substitution model shifts with the K-Pg boundary Phylogenetic logistic regression models show a rapid increase in the probability of a molecular model shift as the time distance to the boundary decreases (Fig. S1a-d).In all models (except those representing mtDNA), time distance to the boundary is a significant predictor of model shift probability (p << 0.05).Considering low (Fig. S1b), mean (Fig. S1c), or high (Fig. S1d) values of phylogenomic discordance has a limited effect on the predicted probability of model shifts.Considering the predicted probability of model shifts across the "merged" nuclear data signal at the time of the oldest shifts in the dataset (marked with leader lines; black curves in Fig. S1), we find a weak inverse relationship between low, mean, and high values of discordance, and model shift probability [63%, 59%, 55%].This implies a very weak confounding effect, with higher levels of discordance leading to a slightly weaker probability of substitution model shifts.These trends are shared with patterns discordance and model shifts in exon and UTR datasets [42%, 32%, 24%, and 44%, 37%, 30%, respectively].The signal of shifts in introns, however, shows the opposite trend, with patterns of low, mean, and high discordance leading to a linear increase in shift probability [19%, 26%, and 34%, respectively], implying that the dynamics of model shifts in introns may be different from those in other nuclear data types.In any case, the bootstrapped confidence intervals for different levels of discordance across introns largely overlap, challenging interpretations.
In sum, and in general, considering phylogenomic discordance in these models does not reject an inference of an association between time distance to the K-Pg boundary and an increased probability of substitution model shifts.Phylogenomic discordance itself is not a significant covariate in any individual pGLM model, and models that exclude discordance entirely are weakly to moderately preferred under AICc (~ 1.5 AICc unit for the "aggregate" model, < 1 AICc units for exons, ~5 AICc units for introns, and ~1.1 AICc units for UTRs).Notably, while there is some collinearity between time distance to the boundary and discordance, the association is only strong in non-phylogenetic regression models (not shown but see Supplementary R code).Thus, models that account for phylogeny may also largely account for collinearity between discordance and time distance to the boundary.

Fig. S1. Model shift probability as a function of time distance to the K-Pg boundary.
Binomial probability of novel molecular model regimes as a function of time distance to the K-Pg boundary (e.g., time distance approaches zero, going from right to left, log scale).Logistic regression curves and associated 70% highest density intervals from bootstrap replicates are depicted over the jittered empirical probability distribution for individual nodes (circular markers).Model projections outside the range of the observed data are shown as dashed curves.In panel A, the probability of a model shift in the aggregate signal from nuclear genetic data (black curve) increases to ~54% at the time of the oldest model shift.The probabilities of shifts across each nuclear genetic data type are somewhat lower (~26%-37%).In contrast, the probability of a shift in mtDNA (all data, or proteins only) is not associated with time distance to the boundary (flat curves).This pattern for mtDNA is consistent with the overall low number of shifts inferred across mtDNA genomes (n < 3).Panels B-D depict logistic regression curves from models that include phylogenomic discordance as a covariate and show patterns of phylogenomic discordance with a white-to-black color ramp.Discordance is not detected as a significant covariate in any model, and models with discordance have higher AICc values, suggesting they are not a better fit to the data.Nevertheless, comparing model predictions with low, mean, and high discordance values shows a very weak confounding effect, with higher levels of discordance leading to a weaker predicted probability of model shifts, except for introns (see supplemental text).

Functional dimensions of sequence variation
We found no statistical evidence that patterns of nucleotide diversity are different across molecular model shift regimes after accounting for non-independence due to phylogeny (F = 0.933 p = 0.98), though we emphasize that this is an exploratory analysis in the absence of population-level sampling at the species level.Next, we found that patterns of synonymous codon usage may be significantly different across model regimes in the omnibus ANOVA, and for many cases of post-hoc comparison after correction for multiple hypothesis tests (SCUO; F: 17.6, p = 0.005, N ̂c; F: 17.23, p = 0.005, N ̂′c F: 27.23, p = 0.0003).These results are depicted in Fig. S2.Notably, the N ̂′c metric appears to be the most sensitive to patterns of variation in codon usage, identifying many cases in which post-hoc pairwise comparisons are statistically different.For example, with N ̂′c the unnamed sister clade to Tinamiformes (a clade uniting Rheiformes, Casuariiformes, and Apterygiformes) appears to have the strongest signal of differentiation, with 5/7 comparisons implying significantly different patterns of codon usage (Fig. S2), sans the groups Passeri or Tinamiformes).In sum, for all eight groups identified by Janus as occupying distinct molecular model regimes in exons, we find evidence of significant differences in patterns of codon usage in comparison to at least one other candidate group.Although it was not our focus here, these results strongly imply that identifying compositional shifts with methods like Janus may be fruitful in generating hypotheses about functional variation in coding sequences.
Interestingly, we find that, after accounting for phylogeny, variation in empirical GC content across all codon positions significantly contributes to taxon partitions identified by Janus.Considering the simple formula: Percentage of Variance Explained by x = (Sum Sq for x / Total Sum of Squares) * 100, we find variation in codon position 1 (47%) contributing similarly to variation in codon position 3 (52%).Codon position 2 contributes a lesser, but not insignificant amount (28%).For positions 1 and 3, the omnibus test reveals a statistically significant result after accounting for phylogeny (p < 0.001).In contrast, the lower amount of explainable variance is reflected in a higher p-value for codon position 2 (p = 0.07).In sum, an examination of the empirical data was consistent with our previous understanding of sequence variation as well as an important role for all codon positions within the coding sequence in the present context.Notably, these patterns (see supplementary R code) are not consistent with an outsized role of the third codon position in driving our results -however, they are consistent with patterns of codon usage as an additional, linked dimension of variation.Assessing the coincidence of molecular shifts and shifts in θ(t) As described in the materials and methods, we assessed statistical support for the hypothesis that shifts in θ(t) may have co-occurred with molecular model shifts using a bootstrapping approach implemented in ℓ1ou (28).Under models relying on the AICc criterion, median support was 87.2% [0.59-0.91]for shifts across the eight-dimensional dataset (pie charts in Fig. 1, Fig. S2).Considering the more conservative pBIC criterion, median support was modestly lower: 76.2% [0- 97.3].For the latter case, three candidate edges received no support relative to the false positive rate: the edge leading to Oscine Passeriformes and the two edges reflecting the stem of Otidae or its unnamed sister clade.Notably, these latter two edges reflect the two deepest splits in Passerea (36), while the edge leading to the MRCA of Passerea receives the strongest relative support under pBIC (97.3%).Applying Fisher's exact test to these inferences supports all cases except those noted above, for which a test of equality in proportion is not rejected at p = 0.05.All other cases are supported at p < 10 -4 under both AICc and pBIC except for two cases: Notopalaeognathae (AICc: p = 0.238, OR = 1.49; pBIC: p = 0.04, OR = 2.1), and the sister clade to Tinamiformes (AICc: p = 0.015, OR = 2.19, pBIC: p < 10 -4 , OR = 4.7) (Fig. S4a, Table S2).
Additionally, unconstrained analyses with ℓ1ou detected 5/12 perfect matches (asterisk inside a circle, Fig. S4b), compared to shifts identified by Janus.All shifts identified by Janus are ~ 1-2 edges away from a shift identified by ℓ1ou in at least one analysis (e.g., pBIC).Thus, we show a tight congruence between patterns identified in molecular data and those identified in life-history data, even when analysis of the latter is not informed by molecular data.We find additional shifts when using ℓ1ou in an unconstrained fashion.However, these shifts are mostly restricted to a temporal interval closely associated with the K-Pg boundary (density plots in S4b).Analyses using the pBIC criterion-the criterion recommended by the authors of ℓ1ou (28)-show the closest correspondence with Janus.The other shifts identified with ℓ1ou are worthy of further examination-however, we omit further discussion of them here as they are of limited relevance in the context of the present manuscript.Fig. S4a.Relative support for shifts in trait optima θ(t).For each criterion (AICc or pBIC), we depict the frequency of empirical positive shift detections from ℓ1ou bootstrapping (light gray) relative to the frequency of false positive detections (dark gray) generated with a dataset simulated under multivariate Brownian motion (e.g., no shifts in θ).Cases in which the empirical positive detection rate is significantly greater than the false positive rate (Fisher's exact test) at p = 0.05 are marked with an asterisk (see materials and methods and Table S2).

Path Density
Edges from Janus to l1ou 0.0 0.2 0.4 Fig. S4b.Unconstrained ℓ1ou analyses.This figure shows results evaluating ℓ1ou and our eightdimensional life-history dataset across three information criteria (AICc, BIC, and pBIC).In the top row, circles mark edges identified across all analyses using Janus, indicating candidate molecular model shift points.Asterisks and phylogeny edge colors indicate shift points identified by ℓ1ou, representing shifts in the multivariate phenotypic optimum θ(t).The second row of the figure shows the temporal sequence of ℓ1ou shifts, indicated by relative density.These plots show that few shifts are detected far from the K-Pg boundary in unconstrained analyses in ℓ1ou.The third row summarizes how many edges separate a given shift point Janus identified to the closest shift identified by ℓ1ou.In addition to 5/12 perfect matches (asterisk inside a circle), all shifts identified by Janus are < 1-2 edges away from a shift identified by ℓ1ou in at least one analysis (e.g., pBIC).

Supplemental models of life-history traits
In the main text, we show visualizations of OUM models in which trait optima θ(t) can shift across specified regimes (Fig. 2, right panel); in this case, macroevolutionary regimes are defined a priori at molecular model shift points.These results stem from an analysis of nuclear genetic data, emphasizing molecular changes associated with the radiation of Neoaves.Indeed, for most cases within this superordinal clade, we observe decreases in the estimates for θ(t) associated with molecular model shifts relative to the ancestral regime θanc.Based on our random forests classifier, we visualized two features with the highest relative variable importance scores ("ChickPC1'' and "adult body mass," Fig. 2, Right).Notably, this set of taxon partitions from Janus includes a plesiomorphic group, uniting our single lineage of Struthio with the clade Galloanserae on the same molecular regime (e.g., Supplemental Fig. 7a).To explore these patterns further, we fit additional OUM models, placing Struthio on a separate regime along with the tree's root.This configuration reflects the aggregate signal of taxon partitions identified across nuclear and mitochondrial datasets (Fig. 1, e.g., adding the shift detected in mitochondrial data for Neognathae).By placing the large-bodied Struthio lineage on a separate regime, we can estimate how its closest sister clades may vary in θ(t) relative to the ancestral regime without being potentially confounded by the inclusion of Galloanseres.For these analyses, we re-estimated OUM models for all numerical lifehistory traits (Fig. S6, also see Fig. S5).
Lastly, Uyeda et al. (151) describes how the analysis of principal components can, under some circumstances, lead to misleading phylogenetic inferences.Specifically, Uyeda et al. (151) describes a phenomenon by which ordination partitions the variance of a multivariate trait such that the major component axes are biased toward "early burst" patterns.As one of the features we examine is a principal component of multivariate developmental mode (40), we checked for this bias by fitting univariate EB, OU, and BM models to each of the eight examined life-history traits and then comparing model fit with AIC weights.We performed these checks using the fitContinuous function in (117) and further validated the results with OUwie (27).
Placing Struthio+root on a separate macroevolutionary regime revealed more complex dynamics within the Palaeognathae (Fig. S6).The Struthio+root regime is detected to have a much larger θ(t), magnifying the apparent shift in life-history trait θ(t) for all other taxon partitions identified by molecular model shift analysis.In this configuration, all candidate regimes identified by Janus show reduced θ(t) for ChickPC1 or adult body mass relative to the ancestral regime.We observe similar patterns for clutch size, and a somewhat weaker signal for generation length (Fig. S5, S6).These analyses corroborate patterns revealed by random forest variable importance analysis (materials and methods), which identified "ChickPC1" and adult body mass as the features that best predict molecular substitution model shifts.Given the challenges in accurately estimating evolutionary parameters like θ(t) from small sample sizes, we present these additional analyses as supporting text.Nonetheless, all estimates are consistent with patterns described in the main text.Further, checks for a bias toward early burst patterns in the data did not reveal any life-history traits for which "EB" was selected as the preferred model."OU" models carried nearly 100% of the model weight in all cases (OU > BM > EB in all cases; see Supplementary R code).Thus, at the scale of the examined phylogeny, we do not see evidence that ordination (in the case of developmental mode data, 40) has biased these data toward an early burst pattern.

Fig. S5.
Complete set of θ(t) shift models (nuclear data configuration).We estimated shifts in θ(t) under an OUM model in OUwie (27).We describe the results for ChickPC1 and adult body mass in the main text, as these two features are identified as being most important with a random forest classifier (See main text and supplementary text).Here, we show the complete set of models for each of the eight continuous life-history traits.Colored distributions reflect parametric bootstrap estimates of θ(t), while grayed-out background distributions reflect a summary of 10,000 simulations of tip values from the model (e.g., time = 0).Notably, analysis of clutch size reveals similar patterns to those observed for ChickPC1 and adult body mass (namely, a reduction in the estimated values for θ(t) for the ancestral regime for many Neoavian clades).Although some fluctuations are observed for other traits, none are as pronounced as for ChickPC1 and adult body mass.

Fig. S6.
Complete set of θ(t) shift models (nuc+mt data configuration).We estimated shifts in θ(t) under an OUM model in OUwie (27).We present ChickPC1 and adult body mass results in the main text, as these two features are identified as most important in a random forest classifier (See main text and supplementary text).Colored distributions reflect parametric bootstrap estimates of θ(t), while grayed-out background distributions reflect a summary of 10,000 simulations of tip values from the model (e.g., time = 0).Here, we show the complete set of models for each of eight continuous life-history traits-with an essential difference from the results presented in Fig. 2 and Fig. S5.Here, an additional partition is included, reflecting a unique molecular model shift identified in the analysis of mtDNA (Neognathae).This has the effect of placing Struthio on a regime separated from Galloanseres (e.g., Fig. 2, Fig. S5).The Galloanseres are thus united with a regime including the Neognathae MRCA.Shift patterns for ChickPC1 and clutch size θ(t) are essentially identical to those presented in previous analyses -however, investigation of body mass reveals more complex dynamics within the Palaeognathae.The Struthio regime contains the tree's root and is naturally estimated to have a much larger θ(t).This, in turn, induces an apparent decrease in θ(t) for body mass in Tinamiformes and more substantial decreases in θ(t) for all other cases.As noted in the supplementary results, these patterns suggest that all molecular model shifts (either nuclear or mtDNA) indicate derived decreases in body mass or ChickPC1 θ(t).Similar fluctuations are observed for other traits (e.g., generation length), but none are as pronounced as the patterns for ChickPC1 and adult body mass.Fig. S7a-d.Maps of estimated equilibrium base frequencies for each taxon partition identified by Janus in analyses of different genetic data types.For each plot a-d, we show a radar plot depicting the relative estimated equilibrium base frequencies for each taxon partition.Colors match those indicated in the legend and correspond to distinct phylogenetic regimes.As noted in the main text and references to Fig. 1, the most significant deviations from empirical base frequencies are identified in our exon dataset (see main text discussion) and are concentrated around a brief interval near the end-Cretaceous mass extinction.All molecular model shifts, except for a shift in exon data in Oscine Passeriformes and mtDNA for Neognathae, are plausibly associated with lineage diversification near the K-Pg boundary (e.g., Fig. S1).

Fig. S8. GC content is strongly linked to life-history variation.
In the above figure (showing linear models), we depict the relationships among average body mass within a particular macroevolutionary regime (labeled markers) identified by Janus against the empirical (left) or estimated equilibrium GC (right) content.We find the expected negative relationship consistent with a GC-biased gene conversion hypothesis in all cases, assuming Ne is broadly and negatively correlated with body mass (materials and methods).In ¾ of the examined cases, most of the observable variance in substitution model parameters (GC%) is explained by its association with life history (R 2 > 0.5).As noted in the main text, the greater deviations between empirical and estimated equilibrium base content observed for exons could be influenced by model fit related to functional constraints, codon usage bias, recombination, or selection (see supplementary analysis of codon usage, Fig. S2, S7a-d).

Statistical performance of Janus
Overall, we observe encouraging performance metrics consistent with those reported by Smith et al. (30).We detected only one false positive out of 1,200 datasets simulated without shifts (Fig. S9a); thus, the rate of false positives was essentially zero across deep or shallow branches.
For datasets simulated with one to four phylogenetically independent shifts (Fig. S9b-e), the average false negative rate declines to negligible levels as dataset size increases and is only significantly greater than zero in a few cases for the smallest simulated datasets of 2kbp (One-Sample t-and Z-tests, n=100, p-values below a significance level of 0.05 after correcting for the false discovery rate ( 102)).For any larger dataset (> 2kpb), the average false negative rate was never significantly larger than zero, though we note that the average false negative rate increases slowly as more shifts are added.For false positives, we observed no cases where the average false positive rate was significantly greater than zero; however, the average false positive rate (and associated standard errors) increases slightly for phylograms based on introns, UTRs, and mtDNA (not exons), as dataset size increases.We suspect this is due to weak over-fitting (30), which is mitigated by only considering strongly supported shifts (e.g., model weights ~ 1.0) in our analyses of empirical data.In any case, as noted above, we emphasize that there are no cases where the average false positive rate is significantly greater than zero.
For cases of nested shifts (Fig. S9f), we observe similarly encouraging performance.In all cases, the average false negative rate declines to negligible levels as dataset size increases.We observe no cases where the average false negative rate is significantly greater than zero (One-Sample t-and Z-tests, n=100, at a significance level of 0.05 after correcting for the false discovery rate).For phylograms based on exons and introns, we observe a weak increase in the average false positive rate as the dataset size increases (to about 10%).Nonetheless, we again observe no cases where the average false positive rate is significantly greater than zero.We again suspect this is due to weak overfitting, as noted above, but we are not concerned that this pattern influences our results based on empirical data, as those analyses only report strongly supported shifts.
In sum, analyses of 9,200 simulated datasets following a wide array of nonhomogeneous patterns of molecular evolution indicate that the performance of Janus is suitable for analyses of substitution model shifts in the context of the Avian phylogeny.In general, under-fitting (e.g., failing to detect simulated shifts) becomes negligible as dataset size increases, and overfitting is never significantly greater than zero with datasets of sufficient size.As all the empirical datasets analyzed in the present study exceed these thresholds, these potential sources of model underperformance should have minimal impact.Our simulations indicate the signal of molecular model shifts described in the main text is robust to the pattern of rapid cladogenesis observed for Aves and that molecular model shifts are only detected when empirical patterns of sequence variation are strong.S9a-f.Summary of 9,200 simulated datasets and associated analyses used to investigate the performance of Janus under a wide range of non-homogeneous model shift scenarios.For each shift scenario, we plot the average false positive and negative rates for each reference phylogram, with associated standard error.Below each plot, we show the results from One-Sample t-tests (Z-tests were similar), which evaluate the null hypothesis that a given configuration and simulated dataset have an average false positive or negative rate that is not greater than zero (after controlling for the false discovery rate due to multiple tests).At a 0.05 significance threshold, few tests suggest that simulated datasets have false positive or negative rates greater than zero.A dashed line at 0.1 (10%) is included for reference.We also provide scaled phylograms to illustrate the range of relative and total branch lengths (tree shapes) evaluated.

Fig. S10. Metabolic scaling parameters are generally predicted by clade mass.
As shown in Fig. 3, we analyzed metabolic data under a Bayesian framework to generate posterior estimates of intercept (β0) and slope (βmass) coefficients from an evolutionary allometric regression model.Here, we depict the posterior estimates for slope and intercept across a fixed allometric shift model recapitulating molecular model regimes across all data types.On the left, we show how metabolic intercept (A, β0) and slope (B, βmass) coefficients are related to the median mass of species in each molecular model regime.Taxa with smaller mass tend to have lower slopes and higher intercept terms (see discussion) than taxa with higher mass.The labels for Datasets 1 and 2 refer to slightly different approaches for estimating group mass.Dataset 1 includes all members within each noted higher taxon, whereas Dataset 2 summarizes the masses of only the "visible" terminal members of each taxon (excluding outliers, e.g., paraphyletic groups, see Fig 1 .).On the right (C), we show the posterior densities for allometric slope (βmass) as a function of allometric intercept (β0).The posterior samples in each region are depicted as a 2D-kernel density with contour lines (see supplementary R scripts in the BMR directory).
For each clade on which we detect a molecular model shift, we note the taxonomic authority and which data types are identified to have model shifts (with an 'x').We report the associated metabolic allometric parameter estimates from bayou for the reader's reference: intercept (β0), slope (βmass), their associated 95% highest posterior density intervals (HPD), the posterior probability of an allometric shift occurring on a particular edge (pp), and the estimated power law describing the association between metabolic rate and body mass ( !!×  With ℓ1ou, we used parametric bootstrapping to estimate the frequency of shift detections and the null false positive rate for each candidate edge under simulated multivariate Brownian motion.We use Fisher's test to ask if the proportion of empirical shift detections in the bootstrapped dataset is higher than that of null false positives in the simulated dataset.All candidate edges are supported under AICc or pBIC criteria, but not both in a few cases.
Table S3.Details of mtDNA dataset assembly.

Fig. S2 .
Fig. S2.Synonymous codon usage patterns are associated with molecular model shifts.For exon data, we evaluated whether patterns of synonymous codon usage differed across phylogenetic regimes identified in the analysis of molecular model shifts.In the top row, we depict matrices of Benjamini-Hochberg-adjusted p-values (102) for post-hoc pairwise comparisons from phylogenetic ANOVA.In the bottom row, we depict T-values for each respective test.See the supplementary text for additional details.

Fig. S3 .
Fig. S3.Patterns of molecular model shifts for each evaluated data type.Shifts in branch color indicate model shifts.Fig.1contains this information in aggregate.Model shifts are most common in exon data, followed by introns, UTRS, and mtDNAs.We note that two model shifts are detected in the analysis of the whole mtDNA dataset (Neognathae and Tinamiformes), whereas analysis of mtDNA proteins or rRNAS separately identifies the same model shift (Neognathae only).We interpret this difference to reflect increased power from the combined mtDNA dataset.See Fig.S7a-dfor an expanded view.

Table S2 . Summary of Fisher's exact test for shifts in life-history trait optima associated with molecular model shifts and higher taxa (36, 150, 152-156).
(146)$ ).Notably, the diverse clade Coraciimoraphae is detected to have molecular model shifts in all nuclear genetic data types and maximal posterior probability for a shift in metabolic allometry.All ESS values were > 500, indicating adequate posterior sampling.Plots showing the evolution of Gelman and Rubin's shrink factor(146)as the number of iterations increases are provided as supplemental data (https://github.com/jakeberv/avian_molecular_shifts/tree/main/BMR).